/**
 * @file   Mesh.h
 * @author Wang Heyu <wang@wukong>
 * @date   Tue May  4 21:45:10 2021
 * 
 * @brief Base class for all types, including irregular and regular.
 * 
 * 
 */

#ifndef __CRAZYFISH_MESH__
#define __CRAZYFISH_MESH__

#include "Point.h"
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <vector>
<<<<<<< HEAD
=======

>>>>>>> 6b22401c3f5f42b7de50a6938c1e67c80a79481a

#define TEMPLATE template <size_t DIM, typename T>

TEMPLATE class Mesh;

TEMPLATE class IrregularMesh;

TEMPLATE
Point<DIM, T> mid_point(const Point<DIM, T> &p0, const Point<DIM, T> &p1);

TEMPLATE
class Mesh
{
public:
	virtual size_t get_n_points() const = 0;
	virtual size_t get_n_edges() const = 0;
	virtual size_t get_n_grids() const = 0;
	virtual size_t get_n_dofs() const = 0;
	virtual const Point<DIM, T> &get_point(size_t _idx) const = 0;
	virtual const Geometry &get_edge(size_t _idx) const = 0;
	virtual const Geometry &get_grid(size_t _idx) const = 0;
	virtual std::vector<int> get_all_boundary();
};

class Easymesh : public Mesh<2, double>
{
private:
	std::valarray<Point<2, double>> nodes;
	std::valarray<Geometry> edges;
	std::valarray<Geometry> grids;
	std::valarray<Point<2, double>> dofs;

public:
	void readData(const std::string &filename);
	Easymesh() : Mesh(){};
	Easymesh(const std::string &filename);
	size_t get_n_points() const;
	size_t get_n_edges() const;
	size_t get_n_grids() const;
	size_t get_n_dofs() const;
	const Point<2, double> &get_point(size_t _idx) const;
	const Point<2, double> &get_dof(size_t _idx) const;
	const std::valarray<Point<2, double>> &get_dofs() const;
	const Geometry &get_edge(size_t _idx) const;
	const Geometry &get_grid(size_t _idx) const;
	Point<2, double> &get_point(size_t _idx);
	Point<2, double> &get_dof(size_t _idx);
	std::valarray<Point<2, double>> &get_dofs();
	Geometry &get_edge(size_t _idx);
	Geometry &get_grid(size_t _idx);
	void build_dofs();
};

const Point<2, double> &Easymesh::get_point(size_t _idx) const
{
	return nodes[_idx];
};

Point<2, double> &Easymesh::get_point(size_t _idx)
{
	return nodes[_idx];
};

const Point<2, double> &Easymesh::get_dof(size_t _idx) const
{
	return dofs[_idx];
};

Point<2, double> &Easymesh::get_dof(size_t _idx)
{
	return dofs[_idx];
};

const std::valarray<Point<2, double>> &Easymesh::get_dofs() const
{
	return dofs;
};

std::valarray<Point<2, double>> &Easymesh::get_dofs()
{
	return dofs;
};

const Geometry &Easymesh::get_edge(size_t _idx) const
{
	return edges[_idx];
};

Geometry &Easymesh::get_edge(size_t _idx)
{
	return edges[_idx];
};

const Geometry &Easymesh::get_grid(size_t _idx) const
{
	return grids[_idx];
};

Geometry &Easymesh::get_grid(size_t _idx)
{
	return grids[_idx];
};

size_t Easymesh::get_n_points() const
{
	return nodes.size();
};

size_t Easymesh::get_n_dofs() const
{
	return dofs.size();
};

size_t Easymesh::get_n_edges() const
{
	return edges.size();
};

size_t Easymesh::get_n_grids() const
{
	return grids.size();
};

/// Most contents copy from Robert Li's AFEPack.
void Easymesh::readData(const std::string &filename)
{
	std::cout << "Reading easymesh data file ..." << std::endl;
	std::cout << "\treading node data ..." << std::flush;
	std::ifstream is((filename + ".n").c_str());
	size_t n_node, n_side, n_element;
	char text[64];

	is >> n_node >> n_element >> n_side;
	is.getline(text, 64);
	nodes.resize(n_node);
	for (size_t i = 0; i < n_node; i++)
	{
		size_t ivtx;
		Geometry::bm_t bm;
		is >> ivtx >> nodes[i][0] >> nodes[i][1] >> bm;
		nodes[i].set_index(ivtx);
		nodes[i].set_boundary_mark(bm);
		nodes[i].set_vertex(0, ivtx);
		nodes[i].set_boundary(0, ivtx);
		nodes[i].set_dof(0, ivtx);
	}
	is.close();

	dofs = nodes;

	std::cout << " OK!" << std::endl;
	std::cout << "\treading side data ..." << std::flush;
	is.open((filename + ".s").c_str());
	size_t n_edge;
	is >> n_edge;
	if (n_edge != n_side)
	{
		std::cerr << "in side file: side number error." << std::endl;
		exit(-1);
	}
	edges.resize(n_side);
	for (size_t i = 0; i < n_side; i++)
	{
		size_t iedg, a, b;
		int ea, eb;
		Geometry::bm_t bm;
		edges[i].set_n_vertex(2);
		edges[i].set_n_boundary(2);
		edges[i].set_n_neighbor(2);
		edges[i].set_n_dofs(2);
		is >> iedg >> a >> b >> ea >> eb >> bm;
		edges[i].set_index(iedg);
		edges[i].set_boundary_mark(bm);
		edges[i].set_vertex(0, a);
		edges[i].set_vertex(1, b);
		edges[i].set_boundary(0, a);
		edges[i].set_boundary(1, b);
		edges[i].set_neighbor(0, ea);
		edges[i].set_neighbor(1, eb);
		edges[i].set_dof(0, a);
		edges[i].set_dof(1, b);
	}
	is.close();
	std::cout << " OK!" << std::endl;

	std::cout << "\treading element data ..." << std::flush;
	is.open((filename + ".e").c_str());
	size_t n_ele;
	is >> n_ele;
	if (n_ele != n_element)
	{
		std::cerr << "in element file: element number error." << std::endl;
		exit(-1);
	}
	size_t n_n;
	is >> n_n;
	if (n_n != n_node)
	{
		std::cerr << "in element file: node number error." << std::endl;
		exit(-1);
	}
	size_t n_s;
	is >> n_s;
	if (n_s != n_side)
	{
		std::cerr << "in element file: side number error." << std::endl;
		exit(-1);
	}
	is.getline(text, 64);
	grids.resize(n_element);
	for (size_t i = 0; i < n_element; i++)
	{
		grids[i].set_n_vertex(3);
		grids[i].set_n_boundary(3);
		grids[i].set_n_neighbor(3);
		size_t igrd, v0, v1, v2, b0, b1, b2;
		int n0, n1, n2;
		is >> igrd >> v0 >> v1 >> v2 >> n0 >> n1 >> n2 >> b0 >> b1 >> b2;
		grids[i].set_index(igrd);
		grids[i].set_boundary_mark(0);
		grids[i].set_vertex(0, v0);
		grids[i].set_vertex(1, v1);
		grids[i].set_vertex(2, v2);
		grids[i].set_neighbor(0, n0);
		grids[i].set_neighbor(1, n1);
		grids[i].set_neighbor(2, n2);
		grids[i].set_boundary(0, b0);
		grids[i].set_boundary(1, b1);
		grids[i].set_boundary(2, b2);
		const Point<2, double> &p0 = nodes[v0];
		const Point<2, double> &p1 = nodes[v1];
		const Point<2, double> &p2 = nodes[v2];
		if ((p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]) < 0)
		{
			std::cerr << "vertices of element " << igrd
					  << " is not correctly ordered." << std::endl;
			// j = g.vertex(2);
			// g.vertex(2) = g.vertex(1);
			// g.vertex(1) = j;
			// j = g.boundary(2);
			// g.boundary(2) = g.boundary(1);
			// g.boundary(1) = j;
		}
		grids[i].set_n_dofs(3);
		grids[i].set_dof(0, v0);
		grids[i].set_dof(1, v1);
		grids[i].set_dof(2, v2);
	}
	is.close();

	std::cout << " OK!" << std::endl;
};

Easymesh::Easymesh(const std::string &filename) : Mesh()
{
	readData(filename);
};

void Easymesh::build_dofs()
{
	int dof_index = 0;
	int n_ele = get_n_grids();
	int n_total_dofs = get_n_points() + get_n_edges();
	dofs.resize(n_total_dofs);
	for (int i = 0; i < n_ele; i++)
	{
		int n_local_dofs = 6;
		grids[i].set_n_dofs(n_local_dofs);
		int n_vtx = grids[i].get_n_vertex();
		for (int k = 0; k < n_vtx; k++)
		{
			Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
			vtx.set_dof(0, -1);
		}
		int n_bnd = grids[i].get_n_boundary();
		for (int k = 0; k < n_bnd; k++)
		{
			Geometry &bnd = get_edge(grids[i].get_boundary(k));
			if (bnd.get_n_dofs() == 2)
			{
				int a = bnd.get_dof(0);
				int b = bnd.get_dof(1);
				bnd.set_n_dofs(3);
				bnd.set_dof(0, a);
				bnd.set_dof(1, b);
				bnd.set_dof(2, -1);
			}
		}
	}
	for (int i = 0; i < n_ele; i++)
	{
		int n_local_dofs = grids[i].get_n_dofs();
		int n_vtx = grids[i].get_n_vertex();
		for (int k = 0; k < n_vtx; k++)
		{
			Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
			int global_dof = vtx.get_dof(0);
			if (global_dof == -1)
			{
				vtx.set_dof(0, dof_index);
				grids[i].set_dof(k, dof_index);
				vtx.set_index(dof_index);
				dofs[dof_index] = vtx;
				dof_index++;
			}
			else
				grids[i].set_dof(k, global_dof);
		}
		int n_bnd = grids[i].get_n_boundary();
		for (int k = 0; k < n_bnd; k++)
		{
			Geometry &edge = get_edge(grids[i].get_boundary(k));
			int global_dof = edge.get_dof(2);
			if (global_dof == -1)
			{
				edge.set_dof(2, dof_index);
				grids[i].set_dof(k + n_vtx, dof_index);
				Point<2, double> &vtx0 = get_point(edge.get_dof(0));
				Point<2, double> &vtx1 = get_point(edge.get_dof(1));
				Point<2, double> vtx = mid_point<2, double>(vtx0, vtx1);
				vtx.set_dof(0, dof_index);
				vtx.set_boundary_mark(edge.get_boundary_mark());
				vtx.set_index(dof_index);
				dofs[dof_index] = vtx;
				dof_index++;
			}
			else
				grids[i].set_dof(k + n_vtx, global_dof);
		}
	}
};

class RegularMesh : public Mesh<2, double>
{
protected:
	Point<2, double> lbc; //x0y0
	Point<2, double> ruc; //x1y1
	size_t nx;
	size_t ny;
	size_t n_dof = 3;

public:
	Point<2, double> node;
	Geometry edge;
	Geometry grid;
	Point<2, double> dof;

	void Initial();
	RegularMesh();
	RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny);
	RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc);
	const Point<2, double> &get_lbc() const { return lbc; };
	const Point<2, double> &get_ruc() const { return ruc; };
	void set_lbc(const Point<2, double> &_lbc) { lbc = _lbc; }
	void set_ruc(const Point<2, double> &_ruc) { ruc = _ruc; }
	size_t get_nx() const { return nx; };
	size_t get_ny() const { return ny; };
	double get_hx() const { return (ruc[0] - lbc[0]) / nx; };
	double get_hy() const { return (ruc[0] - lbc[0]) / ny; };
	void set_nx(size_t _nx) { nx = _nx; };
	void set_ny(size_t _ny) { ny = _ny; };
	void set_nxny(size_t _nx, size_t _ny) { nx = _nx, ny = _ny; };
	size_t get_n_points() const;
	size_t get_n_edges() const;
	size_t get_n_grids() const;
	virtual size_t get_n_dofs() const;
	Point<2, double> &get_point(size_t _idx);
	virtual Point<2, double> &get_dof(size_t _idx);
	Geometry &get_edge(size_t _idx);
	virtual Geometry &get_grid(size_t _idx);
	void _get_grid(size_t _idx);
	std::vector<int> get_all_boundary();
};
void RegularMesh::Initial()
{
	edge.set_n_dofs(2);
	edge.set_n_vertex(2);
	edge.set_n_neighbor(2);
	grid.set_n_neighbor(3);
	grid.set_n_vertex(3);
	grid.set_n_dofs(3);
}
RegularMesh::RegularMesh() : Mesh()
{
	lbc[0] = lbc[1] = 0;
	ruc[0] = ruc[1] = 1;
	nx = ny = 2;
	Initial();
};
RegularMesh::RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny) : Mesh(),
																										 lbc(_lbc), ruc(_ruc), nx(_nx), ny(_ny)
{
	Initial();
}

RegularMesh::RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc) : Mesh(),
																					   lbc(_lbc), ruc(_ruc), nx(2), ny(2)
{
	Initial();
}

Point<2, double> &RegularMesh::get_point(size_t _idx)
{
	size_t ix = (_idx) % (nx + 1);
	size_t iy = (_idx) / (nx + 1);
	double hx = get_hx(), hy = get_hy();
	node[0] = ix * hx;
	node[1] = iy * hy;
	node.set_index(_idx);
	node.set_n_dofs(1);
	node.set_n_vertex(1);
	node.set_vertex(0, _idx);
	if (iy == 0)
		node.set_boundary_mark(1);
	else if (ix == nx)
		node.set_boundary_mark(2);
	else if (iy == ny)
		node.set_boundary_mark(3);
	else if (ix == 0)
		node.set_boundary_mark(4);
	else
		node.set_boundary_mark(0);

	return node;
};

Point<2, double> &RegularMesh::get_dof(size_t _idx)
{
	return get_point(_idx);
};

Geometry &RegularMesh::get_edge(size_t _idx)
{
	size_t ix = (_idx) % (3 * nx + 1);
	size_t iy = (_idx) / (3 * nx + 1);
	edge.set_index(_idx);
	if (ix < nx)
	{
		edge.set_neighbor(0, ((iy == 0) ? -1 : (iy - 1) * (2 * nx) + 2 * ix));
		edge.set_neighbor(1, ((iy == ny) ? -1 : iy * 2 * nx + 2 * ix + 1));
		edge.set_dof(0, iy * (nx + 1) + ix);
		edge.set_dof(1, iy * (nx + 1) + ix + 1);
		edge.set_vertex(0, iy * (nx + 1) + ix);
		edge.set_vertex(1, iy * (nx + 1) + ix + 1);
	}
	else
	{
		ix -= nx;
		if (ix % 2 == 0)
		{
			ix /= 2;
			edge.set_neighbor(0, ((ix == 0) ? -1 : iy * (2 * nx) + 2 * ix - 1));
			edge.set_neighbor(1, ((ix > nx) ? -1 : iy * (2 * nx) + 2 * ix));
			edge.set_dof(0, iy * (nx + 1) + ix);
			edge.set_dof(1, (iy + 1) * (nx + 1) + ix);
			edge.set_vertex(0, iy * (nx + 1) + ix);
			edge.set_vertex(1, (iy + 1) * (nx + 1) + ix);
		}
		else
		{
			ix /= 2;
			edge.set_neighbor(0, iy * (2 * nx) + 2 * ix);
			edge.set_neighbor(1, iy * (2 * nx) + 2 * ix + 1);
			edge.set_dof(0, iy * (nx + 1) + ix + 1);
			edge.set_dof(1, (iy + 1) * (nx + 1) + ix);
			edge.set_vertex(0, iy * (nx + 1) + ix + 1);
			edge.set_vertex(1, (iy + 1) * (nx + 1) + ix);
		}
	}
	if (iy == 0)
		edge.set_boundary_mark(1);
	else if (ix == nx)
		edge.set_boundary_mark(2);
	else if (iy == ny)
		edge.set_boundary_mark(3);
	else if (ix == 0)
		edge.set_boundary_mark(4);
	else
		edge.set_boundary_mark(0);
	return edge;
};

Geometry &RegularMesh::get_grid(size_t _idx)
{
	_get_grid(_idx);
	return grid;
}

void RegularMesh::_get_grid(size_t _idx)
{
	size_t dof[3];
	size_t en[3];
	if (_idx % 2 == 0)
	{
		size_t ix = (_idx / 2) % nx;
		size_t iy = (_idx / 2) / nx;

		dof[0] = iy * (nx + 1) + ix;
		dof[1] = iy * (nx + 1) + ix + 1;
		dof[2] = (iy + 1) * (nx + 1) + ix;

		en[0] = _idx + 1;
		en[1] = (iy == 0) ? -1 : _idx - 1;
		en[2] = (ix == 0) ? -1 : _idx - 2 * nx + 1;
	}
	else
	{
		size_t ix = ((_idx - 1) / 2) % nx;
		size_t iy = ((_idx - 1) / 2) / nx;

		dof[0] = iy * (nx + 1) + ix + 1;
		dof[1] = (iy + 1) * (nx + 1) + ix + 1;
		dof[2] = (iy + 1) * (nx + 1) + ix;

		en[0] = (iy == ny) ? -1 : _idx + 2 * nx - 1;
		en[1] = _idx - 1;
		en[2] = (ix == nx) ? -1 : _idx + 1;
	}
	grid.set_dof(0, dof[0]);
	grid.set_dof(1, dof[1]);
	grid.set_dof(2, dof[2]);

	grid.set_vertex(0, dof[0]);
	grid.set_vertex(1, dof[1]);
	grid.set_vertex(2, dof[2]);

	grid.set_neighbor(0, en[0]);
	grid.set_neighbor(1, en[1]);
	grid.set_neighbor(2, en[2]);

	grid.set_index(_idx);
};

size_t RegularMesh::get_n_points() const
{
	return (nx + 1) * (ny + 1);
};

size_t RegularMesh::get_n_dofs() const
{
	return (nx + 1) * (ny + 1);
};

size_t RegularMesh::get_n_edges() const
{
	return (nx + 1) * ny + (ny + 1) * nx + nx * ny;
};

size_t RegularMesh::get_n_grids() const
{
	return nx * ny * 2;
};

std::vector<int> RegularMesh::get_all_boundary()
{
	std::vector<int> _BndIndex;
	for(int i = 0;i <= nx;i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs()- i - 1);
    }
    for(int j = 1;j <= ny - 1;j++)
    {
        _BndIndex.push_back(j*(nx + 1));
        _BndIndex.push_back((j + 1)*(ny + 1) - 1);
    }
	return _BndIndex;
};

class RegularMeshP2 : public RegularMesh
{
public:
	RegularMeshP2() : RegularMesh() { grid.set_n_dofs(6); };
	RegularMeshP2(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny) : RegularMesh(_lbc, _ruc, _nx, _ny) { grid.set_n_dofs(6); };
	RegularMeshP2(const Point<2, double> &_lbc, const Point<2, double> &_ruc) : RegularMesh(_lbc, _ruc) { grid.set_n_dofs(6); };
	Point<2, double> &get_dof(size_t _idx);
	Geometry &get_grid(size_t _idx);
	size_t get_n_dofs() const;
	std::vector<int> get_all_boundary();
};

Point<2, double> &RegularMeshP2::get_dof(size_t _idx)
{
	double hx = get_hx(), hy = get_hy();
	size_t ix = (_idx) % (2 * nx + 1);
	size_t iy = (_idx) / (2 * nx + 1);
	node[0] = ix * hx / 2;
	node[1] = iy * hy / 2;
	node.set_index(_idx);
	node.set_n_dofs(1);
	node.set_n_vertex(1);
	node.set_vertex(0, _idx);
	if (iy == 0)
		node.set_boundary_mark(1);
	else if (ix == nx)
		node.set_boundary_mark(2);
	else if (iy == ny)
		node.set_boundary_mark(3);
	else if (ix == 0)
		node.set_boundary_mark(4);
	else
		node.set_boundary_mark(0);
	return node;
}

Geometry &RegularMeshP2::get_grid(size_t _idx)
{
	size_t dof[6];
	_get_grid(_idx);
	if (_idx % 2 == 0)
	{
		size_t ix = (_idx / 2) % nx;
		size_t iy = (_idx / 2) / nx;

		dof[0] = 2 * iy * (2 * nx + 1) + 2 * ix;
		dof[1] = 2 * iy * (2 * nx + 1) + 2 * (ix + 1);
		dof[2] = 2 * (iy + 1) * (2 * nx + 1) + 2 * ix;
		dof[3] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 0.5);
		dof[4] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * ix;
		dof[5] = 2 * iy * (2 * nx + 1) + 2 * (ix + 0.5);
	}
	else
	{
		size_t ix = ((_idx - 1) / 2) % nx;
		size_t iy = ((_idx - 1) / 2) / nx;

		dof[0] = 2 * iy * (2 * nx + 1) + 2 * (ix + 1);
		dof[1] = 2 * (iy + 1) * (2 * nx + 1) + 2 * (ix + 1);
		dof[2] = 2 * (iy + 1) * (2 * nx + 1) + 2 * ix;
		dof[3] = 2 * (iy + 1) * (2 * nx + 1) + 2 * (ix + 0.5);
		dof[4] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 0.5);
		dof[5] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 1);
	}
	for (int i = 0; i < 6; i++)
		grid.set_dof(i, dof[i]);
	for (int i = 0; i < 3; i++)
		grid.set_vertex(i, dof[i]);
	return grid;
}

size_t RegularMeshP2::get_n_dofs() const
{
	return (2 * nx + 1) * (2 * ny + 1);
}

std::vector<int> RegularMeshP2::get_all_boundary()
{
	std::vector<int> _BndIndex;
	for(int i = 0;i <= 2 * nx;i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * ny - 1;j++)
    {
        _BndIndex.push_back(j* (2  * nx + 1));
        _BndIndex.push_back((j + 1)*(2 * nx + 1) - 1);
    }
	return _BndIndex;
}

template <typename T>
class RectangleMesh : public Mesh<2, T>
{
private:
	std::valarray<Point<2, T>> nodes;
	std::valarray<Geometry> edges;
	std::valarray<Geometry> grids;
	std::valarray<Point<2, T>> dofs;

public:
	RectangleMesh() : Mesh<2, T>(){};
	RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny);
	virtual size_t get_n_points() const;
	virtual size_t get_n_edges() const;
	virtual size_t get_n_grids() const;
	virtual size_t get_n_dofs() const;
	virtual const Point<2, T> &get_point(size_t _idx) const;
	virtual const Point<2, T> &get_dof(size_t _idx) const;
	virtual const std::valarray<Point<2, T>> &get_dofs() const;
	virtual const Geometry &get_edge(size_t _idx) const;
	virtual const Geometry &get_grid(size_t _idx) const;
};

template <typename T>
const Point<2, T> &RectangleMesh<T>::get_point(size_t _idx) const
{
	return nodes[_idx];
};

template <typename T>
const Point<2, T> &RectangleMesh<T>::get_dof(size_t _idx) const
{
	return dofs[_idx];
};

template <typename T>
const std::valarray<Point<2, T>> &RectangleMesh<T>::get_dofs() const
{
	return dofs;
};

template <typename T>
const Geometry &RectangleMesh<T>::get_edge(size_t _idx) const
{
	return edges[_idx];
};

template <typename T>
const Geometry &RectangleMesh<T>::get_grid(size_t _idx) const
{
	return grids[_idx];
};

template <typename T>
size_t RectangleMesh<T>::get_n_points() const
{
	return nodes.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_dofs() const
{
	return dofs.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_edges() const
{
	return edges.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_grids() const
{
	return grids.size();
};

template <typename T>
RectangleMesh<T>::RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny)
{
	T hx = (_ruc[0] - _lbc[0]) / _nx;
	T hy = (_ruc[1] - _lbc[1]) / _ny;
	nodes.resize((_nx + 1) * (_ny + 1));
	edges.resize(_nx * (_ny + 1) + (_nx + 1) * _ny + _nx * _ny);
	grids.resize(_nx * _ny * 2);

	for (int ix = 0; ix < _nx + 1; ix++)
		for (int iy = 0; iy < _ny + 1; iy++)
		{
			int node_idx = iy * (_nx + 1) + ix;
			nodes[node_idx].set_index(node_idx);
			/// debug
			if (ix == 0 || ix == _nx || iy == 0 || iy == _ny)
				nodes[node_idx].set_boundary_mark(1);
			/// end
			nodes[node_idx].set_vertex(0, node_idx);
			nodes[node_idx].set_boundary(0, node_idx);
			nodes[node_idx].set_dof(0, node_idx);
			nodes[node_idx][0] = ix * hx;
			nodes[node_idx][1] = iy * hy;
		}
	for (int ix = 0; ix < _nx; ix++)
		for (int iy = 0; iy < _ny; iy++)
		{
			int even_ele_idx = (iy * _nx + ix) * 2;

			grids[even_ele_idx].set_n_vertex(3);
			grids[even_ele_idx].set_n_boundary(3);
			grids[even_ele_idx].set_n_neighbor(3);

			grids[even_ele_idx].set_index(even_ele_idx);
			int v0 = iy * (_nx + 1) + ix;
			int v1 = iy * (_nx + 1) + ix + 1;
			int v2 = (iy + 1) * (_nx + 1) + ix;
			grids[even_ele_idx].set_vertex(0, v0);
			grids[even_ele_idx].set_vertex(1, v1);
			grids[even_ele_idx].set_vertex(2, v2);
			grids[even_ele_idx].set_n_dofs(3);
			grids[even_ele_idx].set_dof(0, v0);
			grids[even_ele_idx].set_dof(1, v1);
			grids[even_ele_idx].set_dof(2, v2);

			int odd_ele_idx = (iy * _nx + ix) * 2 + 1;

			grids[odd_ele_idx].set_n_vertex(3);
			grids[odd_ele_idx].set_n_boundary(3);
			grids[odd_ele_idx].set_n_neighbor(3);

			grids[odd_ele_idx].set_index(odd_ele_idx);
			v0 = (iy + 1) * (_nx + 1) + ix + 1;
			v1 = (iy + 1) * (_nx + 1) + ix;
			v2 = iy * (_nx + 1) + ix + 1;
			grids[odd_ele_idx].set_vertex(0, v0);
			grids[odd_ele_idx].set_vertex(1, v1);
			grids[odd_ele_idx].set_vertex(2, v2);
			grids[odd_ele_idx].set_n_dofs(3);
			grids[odd_ele_idx].set_dof(0, v0);
			grids[odd_ele_idx].set_dof(1, v1);
			grids[odd_ele_idx].set_dof(2, v2);
		}
	dofs = nodes;
};

TEMPLATE
Point<DIM, T> mid_point(const Point<DIM, T> &p0, const Point<DIM, T> &p1)
{
	Point<DIM, T> mid_pnt;
	for (int i = 0; i < DIM; i++)
		mid_pnt[i] = (p0[i] + p1[i]) * 0.5;
	return mid_pnt;
};

#undef TEMPLATE

#else
// DO NOTHING.
#endif
